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Abstract 


Improved  factors  of  safety  for  quantitative  estimates  for  grid  and  time 
convergence  uncertainties  for  CFD  solutions  are  proposed  for  situations  when  Richardson 
extrapolation  estimated  order  of  accuracy  pk  is  larger  than  the  theoretical  order  of 

accuracy  pk  and  correction  factor  l<C/t<2.  The  improved  uncertainty  estimates  are 

shown  to  provide  more  reasonable  intervals  of  uncertainty  for  pk  >  pk  (l<C/(<2). 


I.  INTRODUCTION 


Current  procedures  for  quantitative  estimates  for  grid  and  time  convergence 
uncertainties  for  CFD  solutions  are  based  on  first-order  Richardson  extrapolation  (RE) 
error  estimates  with  factors  of  safety  ( FS)  used  for  expanded  uncertainty  estimates. 
Roache  [1,2]  proposed  the  grid-convergence  index  (GCI)  with  FS=\ .25  for  systematic 
grid-triplet  studies  using  RE  estimate  for  order  of  accuracy  pk  ( k=G  for  grid,  k=T  for 
time,  and  k=P  for  similar  parameters)  and  FS=  3  for  2-grid  sensitivity  studies  using 
theoretical  estimate  for  order  of  accuracy  pk  .  The  GCI  is  widely  used  and  recommended 

by  ASME  [3]  and  AIAA  [4].  The  authors  and  colleagues  proposed  a  correction  factor 
(Ca)  method  [5,6]  with  linearly  increasing  FS  vs.  distance  from  the  asymptotic  range 

(AR)  Ck  =  (rkk  -1)1  (rkM  - 1) ,  which  was  based  on  analytical  benchmarks  that  approach 

the  AR  with  pk  <  pk  .  FS  was  reflected  for  pk  >  pk  .  FS  using  Ck  is  smaller  than 

FS=\  .25  for  solutions  very  close  to  the  AR,  whereas  FS  using  Ck  is  much  larger  than 
FS=  \  .25  for  solutions  far  from  the  AR,  which  is  the  typical  situation  for  industrial 
applications.  Ck  has  the  “common-sense”  advantage  compared  to  GCI  in  providing  a 
quantitative  metric  to  determine  proximity  of  the  solutions  to  the  AR  and  approximately 
accounts  for  the  effects  of  higher-order  RE  tenns.  The  Ck  method  has  been  used  in  ship 
hydrodynamics  CFD  workshops. 

A  deficiency  of  both  GCI  and  Ck  methods  is  that  for  pk  >  p,  the  uncertainty 
estimates  are  unreasonably  small  in  comparison  to  uncertainty  estimates  for  pk  <  pk 
with  similar  grid  refinement  ratio  rk ,  solution  changes  between  fine  and  medium 
grid/time  ek  and  distance  \pk  —  pm |  from  the  AR,  which  results  from  too  small  RE  error 
estimate: 


(1) 
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In  the  GCI  method  FS=  \  .25  is  much  too  small  and  in  the  Ca  method  linearly  increasing 
FS  is  also  too  small.  Herein,  an  improvement  to  the  Ca  method  is  proposed  with 
polynomial  increasing  FS  for  pk  >  pk  based  on  reflecting  the  grid/time  uncertainty  from 

pk  <  pk  for  pk  >  pk ;  .  The  improved  uncertainty  estimates  are  shown  to  provide  more 

reasonable  intervals  of  uncertainty  for  pk  >  pk  . 


II.  IMPROVED  FS  FOR  RE  ESTIMATED  LARGER  THAN  THEORETICAL 

ORDER  OF  ACCURACY 


Grid  and  time  convergence  studies  are  conducted  with  multiple  solutions  (at  least 
3)  using  systematically  refined  grid  sizes  or  time  steps.  For  monotonic  convergence, 
procedures  for  estimating  grid  and  time  errors  are  based  on  RE,  which  assumes  that  the 
error  terms  are  in  the  form  of  power  series  expansion.  Results  from  the  numerical 
solution  of  the  one-dimensional  wave  and  two-dimensional  Laplace  equation  analytical 
benchmarks  show  that  RE  error  estimate  using  Eqn.  (1)  has  the  right  form/trends,  but  is 
only  qualitatively  not  quantitatively  accurate  due  to  poorly  estimated  pk  <  pk 


„  _Hede* J 

Pk  ~  ,  t  \ 
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(2) 


except  when  solutions  are  very  close  to  the  AR.  The  error  estimate  can  be  improved  using 

correction  factors,  i.e.,  =  CkS*RE  ,  which  is  used  to  estimate  the  uncertainty  for 
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uncorrected  solutions  by  bounding  the  error  Sk  by  the  sum  of  the  absolute  value  of  the 
corrected  RE  error  estimate  and  the  absolute  value  of  the  amount  of  the  correction  with 


provision  for  10%  FS  in  the  limit  of  Ck  =1.  FS  =  Uk 


'RE,. 


is  reflected  from  pk  <  pk 


for  pk  >  pk  [6].  Thus  for  uncorrected  solutions, 
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0.875  <  Ck  <1.125 
0<CA<  0.875  or  Ck  >1.125 


(3) 


Uk=FS 


[9.6(1-C,)!  +  1 

.1 

^REh 

[2|1-C,|  +  1]|^ 

For  corrected  solutions,  Uk  is  based  on  the  absolute  value  of  the  amount  of  the 
correction: 


uk=\ 


[2.4(1  —  Ckf  +  0.1 


0.75  <C,  <1.25 


0  <  Q  <  0.75  or  C,  >  1 .25 


(4) 


Ck  [6]  method  is  equivalent  to  the  GCI,  but  with  a  variable  FS  that  increases  linearly 
with  the  distance  of  solutions  increases  from  the  AR,  as  shown  in  Fig.  1. 

Verification  studies  have  shown  that  the  estimates  of  U k  and  U k  using  Eqns.  (3) 

and  (4)  are  too  conservative  for  pk  >  pk  ,  as  explained  previously.  An  improved 
approach  is  to  reflect  the  uncertainty  itself  with  respect  to  the  distance  from  the  AR 
instead  of  reflecting  the  FS.  First,  rkk  in  Eqn.  (1)  is  re-expressed  based  on  the  definition 
of  Cp 

<k  =Ck(rkt,h  -])  +  1  (5) 

Second,  Eqn.  (1)  is  substituted  into  Eqn.  (3)  for  0<Ca<0.875  with  the  use  of  Eqn.  (5), 
which  results  in  an  alternative  fonn  of  Uk 


tff=[2(l-C,)  +  l] 


0<Q  <0.875 


(6) 
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Figure  1.  Factors  of  safety  for  correction  factor  and  GCI  verification  methods. 


Third,  to  enforce  the  same  F4  for  1 . 1 25<Q<2  as  the  L4  at  the  same  distance  to  Ca=1 
within  0<Ca<0.875,  Ca  in  Eqn.  (6)  is  replaced  by  2-Ca.  Thus  for  the  same  rk,  pk  ,  and 
ek  ,  Equation  (6)  becomes: 


[2(Ct-l)  +  l] \SL  1.125  <Ck  <2 


Additionally,  two  3rd  order  polynomials  instead  of  two  quadratic  functions  as  used  in  [6] 
are  used  to  generate  smoother  curves  of  FS  for  0.875<Ca<1.0  and  1.0<C*<1.125, 
respectively.  The  use  of  higher  order  polynomials  allows  not  only  the  FS  magnitude  but 
also  the  first  order  derivative  of  FS  with  respect  to  Ca  to  be  continuous  at  Ca=0.875,  1.0, 
and  1.125,  which  are  used  to  determine  the  four  unknown  coefficients  for  each 
polynomial.  Following  [6],  the  magnitude  and  1st  order  derivative  of  FS  at  Ca=1  are  1.1 
and  0,  respectively.  Incorporating  all  of  the  above  revisions  the  Uk  and  Uk  are  given  by: 
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Uk=\ 


[2(l-Q)  +  l] 

[-25.6(1  -  Ckf  +12.8(1  -Ckf  +l.l] 
[-135.8(C*  -l)3  +49.4(C,  -l)2  +  l.l] 


'REk, 


c„ 


2  —  C, 


{2(c,-i)+i]U; 


0  <  C,  <  0.875 


0.875  <C,  <1.0 


1.0 <C,  <1.125 


1.125  <Ct  <2 


(8) 


uk  = 


'REh 


(1  -Ck) 

[-3.2(l-C,)3+3.2(l-C,)2+0.l] 
[—16.98(0*  - 1)3  + 12.35(0*  - 1)2  +  0.  l] 


'REt, 


'RE. 


rCl  +  2Ct-3'' 
3  -C, 


0  <  C  <  0.75 


0.75  <C  <1.0 


1.0 <C„  <1.25 


1.25  <Ck  <2 


(9) 


Compared  to  Ck  [6],  the  improved  Ck  method  introduces  an  additional  term  Ck /  ( 2  -  Ck ) 
to  compute  Uk  for  1 . 1 25<Ck<2.  When  Ck  increases,  this  tenn  increases  rapidly  from  1  to 

infinity,  which  amplifies  FS  when  solutions  are  further  away  from  the  AR.  The  improved 
Ck  method  is  only  applicable  for  0<Ca<2.  C7=0  is  the  border  of  convergence  and 
divergence  such  that  grid  errors/uncertainties  are  infinite  due  to  infinite  <7*^  as  a  result 

of  /?a=0,  i.e.  solution  changes  for  the  medium  and  fine  grids  are  equal  to  those  for  the 
coarse  and  medium  grids.  For  C/(>2,  solutions  are  too  far  from  the  AR  and  also  regarded 
as  divergent.  Figure  1  compares  FS  predicted  by  the  improved  Ck,  Ck  [6],  and  GCI 
methods  with  a  zoomed  in  view  near  the  AR  shown  in  Fig.  2. 

Compared  with  the  Ck  [6]  method,  the  improved  Ck  method  is  more  conservative 
for  C/(>0.875  except  with  the  same  FS  at  C/=  I .  The  intersection  points  between  the 
improved  Ck  and  GCI  methods  depends  on  value  FS  used  in  GCI,  e.g.,  for  FS=  \  .25 
intersection  points  are  Ck=  (0.875,  1.06)  and  (0.75,  1.12)  for  uncorrected  and  corrected 
solutions,  respectively.  When  solutions  are  between  the  intersection  points,  i.e.,  closer  to 
the  AR,  the  GCI  method  is  more  conservative  than  the  improved  Ck  method.  When 
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solutions  are  outside  the  intersection  points,  i.e.,  further  from  the  AR,  the  GCI  method  is 
less  conservative  than  the  improved  Q  method. 


0.9  095  1  1  05  1  1  1.15 


Figure  2.  A  zoomed  in  view  of  factors  of  safety  for  correction  factor  and  GCI 
verification  methods. 

III.  EXAMPLES  FOR  SHIP  HYDRODYNAMICS  APPLICATIONS 

To  demonstrate  that  the  improved  C*  method  predicts  more  reasonable  intervals 
of  grid  uncertainties  than  Ca  [6]  and  GCI  methods  for  industrial  applications,  the  three 
methods  are  applied  for  a  recent  study  [7]  that  used  computational  towing  tank 
procedures  for  single  run  curves  of  resistance  and  propulsion  for  the  high-speed  transom 
ship  Athena  barehull  with  a  skeg  using  the  general-purpose  solver  CFDShip-Iowa-V.4 
[8].  Extensive  verification  (for  grid)  and  validation  (not  shown  herein)  studies  are 
conducted  by  continuously  refining  the  grid  from  the  coarsest  grid  (grid  7  with  360,528 
points;  y[  =4.26)  to  the  finest  grid  (grid  1  with  8.1  million  points;  \\  =1.52)  for  the 
Athena  bare  hull  with  skeg  with  2  degrees  of  freedom  (pitch  and  heave)  at  Froude 
number  (Fr)  0.48.  The  grids  are  designed  with  a  systematic  grid  refinement  ratio  re  =  21/4, 
which  allows  9  sets  of  grids  for  verification  and  validation  (V&V)  with  5  sets  with  re  = 
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21/4  (5,6,7;  4,5,6;  3,4,5;  2,3,4;  and  1,2,3),  3  sets  with  rc  =  2m  (3,5,7;  2,4,6;  and  1,3,5),  and 
1  set  with  i'g  =  23  4  (1,4,7).  Figure  3(a)  and  3(c)  show  the  solutions  with  EFD  data  for  the 
resistance  coefficients  and  ship  motions,  respectively.  Figure  3(b)  and  3(d)  show  the 
relative  solution  changes  between  two  successive  grids  with  iterative  errors  for  the 
resistance  coefficients  and  ship  motions,  respectively.  The  coarsest  grid  7  is  too  coarse  as 
its  solution  is  out  of  the  trend  shown  for  the  other  grid  solutions.  eN  shows  systematic 

decreasing  for  CTx,  Qx,  and  trim  with  UI  <  eN  for  the  coarse  grids.  eN  shows  oscillatory 
decreasing  for  CPX  and  sinkage,  which  is  caused  by  the  problem  of  separating  the  iterative 
errors  Ui  and  eN  for  the  tine  grids  as  they  are  of  the  same  order  of  magnitude.  Overall  Ui 

is  insensitive  to  the  refinement  of  grids  and  the  average  U ,  is  0.18%,  0.15%,  and  0.2% 
for  CTX,  sinkage  and  trim,  respectively.  CJX  monotonically  converges  for  6  sets  of  grids 
except  those  with  the  coarsest  grid  7  involved,  whereas  motions  are  more  difficult  to 
converge.  Verification  results  for  monotonic  converged  solutions  are  presented  in  Tables 
1  and  2  for  the  total  resistance  coefficient  and  ship  motions,  respectively.  C(;  shows  large 

range  of  oscillations  (0.07<CG<2.42  for  CTX  and  0.40<CG<  16.92  for  ship  motions) 
indicating  that  the  solutions  are  not  yet  in  the  AR. 

Table  1.  Verification  study  for  Ctx  of  Athena  bare  hull  with  skeg  (Fr=0.48).  Ug  is  %Sfine; 
Ctx  is  based  on  static  wetted  area;  Factor  of  safety  for  GCI  is  1.25 


Grids 

Refinement 

Ratio 

Rg 

{£G2l  /fG32) 

Pg 

Cc 

UG  (%) 

Xing  and 
Stem 

Ck[  6] 

GCI 

2,4,6 

0.63 

1.32 

0.58 

4.90 

4.90 

3.34 

1,3,5 

V2 

0.40 

2.66 

1.51 

3.59 

1.16 

0.72 

4,5,6 

0.97 

0.16 

0.07 

125.2 

125.2 

52.7 

3,4,5 

yfl 

0.80 

1.27 

0.59 

7.23 

7.23 

4.98 

2,3,4 

yfl 

0.60 

2.98 

1.64 

8.73 

4.27 

1.07 

1,2,3 

0.50 

4.00 

2.42 

— 

1.11 

0.58 

As  shown  in  Table  1,  Ug  of  C'TX  for  grids  (4,5,6)  is  unreasonable  large  as  it  is  too 
far  away  from  the  asymptotic  range.  Ug  of  CJX  for  grids  (1,2,3)  using  the  C*  [6]  is 


7 


unreasonable  small  due  to  the  deficiency  discussed  above.  Excluding  these  two  numbers, 
the  average  Ug  are  2.53%,  4.39%,  and  6.11%  for  the  GCI,  Ca  [6],  and  improved  C* 
methods,  respectively.  For  pG  >  pG  on  grids  (2,3,4),  the  improved  Ca  method  predicts 

more  reasonable  Uq  (8.73%),  which  is  2  times  the  magnitude  using  Ca  [6]  method  and 
one  order  of  magnitude  larger  than  that  using  GCI.  When  solutions  are  closer  to  the  AR, 
the  differences  between  the  Ug  using  the  three  methods  decrease.  As  shown  in  Table  2 
for  trim  on  grids  (1,3,5)  where  Cg=1.09,  U(j  is  of  the  same  order  of  magnitude  for  the 
three  methods.  Nonetheless,  the  improved  Ca  method  is  more  conservative  than  GCI  and 
GCI  is  more  conservative  than  the  Ca  [6]  method,  which  is  consistent  with  the 
observation  in  Fig.  2  for  1.06<Cg^1.125. 


Table  2.  Verification  study  for  motions  of  Athena  bare  hull  with  skeg  (0=0.48).  Ug  is 
%Sfme;  Factor  of  safety  for  GCI  is  1.25 


Parameter 

Grids 

Refinement 

Ratio 

Rg 

(c?,,  /c?32) 

Pg 

CG 

Ua(%) 

Xing  and 
Stern 

Ck[  6] 

GCI 

Sinkage 

1,3,5 

0.31 

3.4 

2.25 

— 

1.80 

0.64 

Sinkage 

2,3,4 

0.13 

12 

16.92 

— 

1.37 

0.05 

Trim 

1,3,5 

0.48 

2.13 

1.09 

4.67 

3.88 

4.12 

Trim 

4,  5,6 

V 2 

0.86 

0.89 

0.40 

42.87 

42.87 

24.42 

Trim 

2,3,4 

^2 

0.53 

3.69 

2.16 

— 

8.91 

3.35 

Trim 

1,2,3 

^2 

0.53 

3.71 

2.18 

— 

4.64 

1.73 

Overall  the  improved  Ca  method  provides  more  reasonable  uncertainty  estimates 
for  pG  >  pGi  than  the  Ca  [6]  and  GCI  methods.  More  accurate  and  efficient  iterative 

methods  (e.g.  multi-grid)  are  needed  to  speed  up  the  convergence  and  reduce  the  Ur 
especially  for  the  fine  grids  for  improved  assessment  of  grid  convergence.  Further 
refinement  with  yf  <1  may  also  help  reach  the  AR  but  requires  at  least  38  million  grid 
points,  which  raises  issues  of  code  efficiency  and  available  computer  resources. 


esistan 


Grids  (N) 


(c)  (d) 

Figure  3.  Verification  for  resistance  and  motions  for  Athena  bare  hull  with  skeg 
(Fr=0.48):  (a)  resistance  coefficients,  (b)  relative  change  =|(.S'V_1  —  ,S'V ) /.S’,  | x  1 00  and 

iterative  errors  for  resistance  coefficients,  (c)  sinkage  and  trim,  (d)  relative  change  c  y  and 
iterative  errors  for  sinkage  and  trim. 

IV.  CONCLUSIONS  AND  FUTURE  WORK 

Improved  factors  of  safety  for  quantitative  estimates  for  grid  and  time  convergence 
uncertainties  for  CFD  solutions  are  proposed  for  situations  when  Richardson 
extrapolation  estimated  order  of  accuracy  pk  is  larger  than  the  theoretical  order  of 
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accuracy  pki  .  The  improved  uncertainty  estimates  are  shown  to  provide  more  reasonable 
intervals  of  uncertainty  for  pk  >  pk  and  l<C/(<2. 

Since  numerical  solutions  of  the  analytical  benchmarks  conducted  so  far  approach 
the  AR  with  pk  <  pk  ,  it  is  desirable  to  validate  current  verification  procedures  using 

more  advanced  numerical  benchmarks  for  complex  flows  such  as  backward-facing  step 
flow  as  the  solutions  will  likely  approach  the  AR  similarly  as  industrial  applications,  i.e., 
with  oscillatory  1  -Q>  The  more  advanced  numerical  benchmarks  can  also  be  used  to 
validate  the  validation  procedures. 
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